suppressMessages(library(Seurat))
suppressMessages(library(ggplot2))
suppressMessages(library(SingleCellExperiment))
sce <- readRDS("/mnt/nmorais-nfs/marta/pB_joana/pC_data/sce-mono-macs-young-without-residents.rds")
sce
## class: SingleCellExperiment 
## dim: 16864 12313 
## metadata(0):
## assays(5): counts logcounts limma limma_lane limma_sortingday
## rownames(16864): Xkr4 Mrpl15 ... CAAA01147332.1 AC149090.1
## rowData names(0):
## colnames(12313): AAACCCAAGGGTTAAT-1 AAACCCACAAGTACCT-1 ...
##   TTTGATCGTTGCTCAA-1 TTTGATCGTTTCTTAC-1
## colData names(18): total_counts log10_total_counts ... doublet
##   cell_type
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

Create new seurat object

norm_mat <- (2**assays(sce)$limma)-1
ln_mat <- log(norm_mat + 1)
srt <- CreateSeuratObject(counts = ln_mat, min.cells = 0, min.features = 0)
## Warning: Non-unique features (rownames) present in the input matrix, making
## unique
## Warning: Non-unique cell names (colnames) present in the input matrix, making
## unique
srt
## An object of class Seurat 
## 16864 features across 12313 samples within 1 assay 
## Active assay: RNA (16864 features, 0 variable features)
srt <- FindVariableFeatures(srt, selection.method = "vst", nfeatures = 3000)
# Identify the 20 most highly variable genes
top20 <- head(VariableFeatures(srt), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(srt)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
plot2
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 49 rows containing missing values (geom_point).

all.genes <- rownames(srt)
srt <- ScaleData(srt, features = all.genes)
## Centering and scaling data matrix
srt <- RunPCA(srt, features = VariableFeatures(object = srt), npcs = 100)
## PC_ 1 
## Positive:  Itm2b, Apoe, C1qc, C1qb, Selenop, Ms4a7, C1qa, Timp2, Pltp, Trem2 
##     Hexa, Hexb, Tmem86a, Rgs10, Arhgap22, Tmcc3, Ctsh, Ctsb, Tnfaip8l2, Frmd4b 
##     Syngr1, Fcrls, Aif1, Serpinb6a, Mertk, Otulinl, Irf8, Itga9, Cxcl16, Maf 
## Negative:  Clec4e, Pim1, Fn1, Clec4d, Emilin2, Cd14, Ahnak, S100a6, Ehd1, Rab11fip1 
##     Cebpb, Txnrd1, S100a4, Gda, F10, Pde4b, Nfkb1, Lilrb4a, Ppp4r2, Thbs1 
##     Cd44, Malt1, Slc7a11, Slfn2, S100a10, Sod2, Adora2b, Plcb1, Ccl9, Srgn 
## PC_ 2 
## Positive:  Ctsd, Fabp5, Ctsl, Cd63, Spp1, Ftl1, Ctsb, Lhfpl2, Pf4, Gpnmb 
##     Lgals1, Vat1, Trem2, Emp1, Cd68, Fth1, Arhgap10, Basp1, Gde1, Cstb 
##     Slc6a8, Syngr1, Adam8, Cd36, Lgals3, Il7r, Cpd, Abca1, Fnip2, Timp2 
## Negative:  Ly6e, Ifi209, Ms4a4c, Samhd1, Plac8, Mndal, Ifi213, A330040F15Rik, Ifitm3, Ifi203 
##     Txnip, Stat2, Ifi204, Irf7, Rnf213, Dleu2, Slfn1, Ccr2, Ifi211, Slfn5 
##     Ifi47, Trim30a, Parp14, Isg15, Phf11b, Ifitm6, Plbd1, Arhgap15, Cd74, Ccnd3 
## PC_ 3 
## Positive:  Cytip, Xylt1, Fyn, Runx3, Napsa, Csf2rb, Plbd1, Lsp1, Msrb1, Gngt2 
##     Itgal, Gan, Cd52, Mir142hg, Hp, Gcsh, Gsr, Antxr2, Gpr141, Pde3b 
##     Runx2, Tspan13, Lmnb1, Itgax, Cbfa2t3, Gpr132, Mcemp1, Arl4c, Actg1, Gcnt2 
## Negative:  Tnf, Cxcl2, Ccl2, Ccl7, Nfkbiz, Il10, Tnfaip3, Cxcl1, Nfkbia, Phlda1 
##     Gm15726, Egr1, Marcksl1, Mrc1, Sdc4, Gadd45b, Sertad2, Mir155hg, Tnfsf9, Sash1 
##     Errfi1, Orai2, G530011O06Rik, Arhgap10, Maff, Gm14636, Mthfs, Pf4, Ccl12, Nrp2 
## PC_ 4 
## Positive:  Zswim6, Pid1, Trps1, Camk1d, Cmss1, Zfp608, Ssh2, Nedd9, Dleu2, Pecam1 
##     Grk5, Gm26917, Foxp1, Ints6, Runx1, Zfp710, Fbxl18, Pde4c, Gm26887, Ccnb1ip1 
##     Slc9a9, Neat1, Sgms1, Dand5, A530013C23Rik, Pdpk1, Ptprj, Ccnd3, Slc12a6, Tmem164 
## Negative:  Isg15, Rsad2, Irf7, Bst2, Phf11b, Ifit3, Oasl2, Ifi211, Zbp1, Ifit2 
##     Phf11d, Slfn5, Ms4a4c, Oas3, Ifi47, Ifi204, Usp18, Ifit1, Lgals3bp, Trim30a 
##     Ifi209, Mndal, Znfx1, A330040F15Rik, Ifi206, Stat1, Parp14, Irgm1, Isg20, Ifi203 
## PC_ 5 
## Positive:  Nfe2l2, H2-Aa, Tmem176b, H2-Ab1, Rel, Il1b, Cd74, Runx3, Ccrl2, Pde4b 
##     Crem, H2-Eb1, Trf, Tmem176a, Jmy, Cd83, Nr4a3, Lpcat2, Abca1, H2-DMb1 
##     Tspan13, Fcgr2b, H2-DMa, Kdm6b, Arl5c, Apoe, Itga9, Gpr132, Cxcl16, Mpp7 
## Negative:  Capg, Anxa2, S100a10, S100a6, S100a4, Lgals1, Anxa1, Crip1, Lpl, Vim 
##     Lgals3, Pkm, Tagln2, Calr, Rnh1, Pdia6, Tmsb10, Msr1, Adam8, Ccl9 
##     S100a11, Qpct, Cdk14, Cd180, Cd36, Cd68, Psd3, Hsp90b1, Vat1, Slfn4
ElbowPlot(srt, ndims = 100) + 
  geom_vline(xintercept = 25, color = "red", alpha = 0.5) + 
  theme(text = element_text(size = 10), 
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10))

srt$sample <- sce$sample
srt$lane <- sce$lane
srt$sorting_day <- sce$sorting_day
srt$cell_type <- sce$cell_type
DimPlot(srt, reduction = "pca", group.by = "cell_type")

srt <- RunUMAP(srt, n.components = 10, features = VariableFeatures(srt))
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 17:06:23 UMAP embedding parameters a = 0.9922 b = 1.112
## 17:06:23 Read 12313 rows and found 3000 numeric columns
## 17:06:23 Using Annoy for neighbor search, n_neighbors = 30
## 17:06:23 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 17:06:33 Writing NN index file to temp file /tmp/RtmpaxeuJs/filed94337e35537
## 17:06:33 Searching Annoy index using 1 thread, search_k = 3000
## 17:08:59 Annoy recall = 100%
## 17:09:00 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 17:09:01 Initializing from normalized Laplacian + noise (using irlba)
## 17:09:02 Commencing optimization for 200 epochs, with 667944 positive edges
## 17:09:12 Optimization finished
srt <- RunTSNE(srt,
             dim.embed = 3,
            dims = 1:25)
saveRDS(srt, "/mnt/nmorais-nfs/marta/pB_joana/pC_data/srt-mono-macs-young-without-residents-bec.rds")
table(srt$cell_type,srt$sample)
##                           
##                             yg0  yg1  yg3  yg5
##   CD11a Myeloid              89    7    9    6
##   Infiltrating Type 1       293 1371  314   46
##   Infiltrating Type 2        26 1865   46   12
##   Macrophages intermediate   15   25 2784 2131
##   Macrophages Repair         15   12 1254 1922
table(srt$sorting_day,srt$sample)
##       
##         yg0  yg1  yg3  yg5
##   day1    0 3280 4407    0
##   day2  438    0    0 4117
table(srt$lane,srt$sample)
##        
##          yg0  yg1  yg3  yg5
##   lane1    0 3280    0 4117
##   lane2  438    0 4407    0
table(srt$cell_type,srt$sorting_day)
##                           
##                            day1 day2
##   CD11a Myeloid              16   95
##   Infiltrating Type 1      1685  339
##   Infiltrating Type 2      1911   38
##   Macrophages intermediate 2809 2146
##   Macrophages Repair       1266 1937
table(srt$cell_type,srt$lane)
##                           
##                            lane1 lane2
##   CD11a Myeloid               13    98
##   Infiltrating Type 1       1417   607
##   Infiltrating Type 2       1877    72
##   Macrophages intermediate  2156  2799
##   Macrophages Repair        1934  1269
colors <- c("darkorange4","red", "yellow", "green4", "blue" )
DimPlot(srt, reduction = "tsne", group.by = "cell_type", cols = colors, pt.size = 0.5)

DimPlot(srt, reduction = "tsne", group.by = "sample", pt.size = 0.5)

DimPlot(srt, reduction = "tsne", group.by = "lane", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt, reduction = "tsne", group.by = "sorting_day", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$lane == "lane1"], reduction = "tsne", group.by = "lane", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$lane == "lane2"], reduction = "tsne", group.by = "lane", cols = c("blue"), pt.size = 0.5)

DimPlot(srt[,srt$sorting_day == "day1"], reduction = "tsne", group.by = "sorting_day", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$sorting_day == "day2"], reduction = "tsne", group.by = "sorting_day", cols = c( "blue"), pt.size = 0.5)

DimPlot(srt, reduction = "umap", group.by = "cell_type", cols = colors, pt.size = 0.5)

DimPlot(srt, reduction = "umap", group.by = "sample", pt.size = 0.5)

DimPlot(srt, reduction = "umap", group.by = "lane", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt, reduction = "umap", group.by = "sorting_day", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$lane == "lane1"], reduction = "umap", group.by = "lane", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$lane == "lane2"], reduction = "umap", group.by = "lane", cols = c("blue"), pt.size = 0.5)

DimPlot(srt[,srt$sorting_day == "day1"], reduction = "umap", group.by = "sorting_day", cols = c("orange", "blue"), pt.size = 0.5)

DimPlot(srt[,srt$sorting_day == "day2"], reduction = "umap", group.by = "sorting_day", cols = c( "blue"), pt.size = 0.5)

DimPlot(srt[,srt$sample == "yg0"], reduction = "umap", group.by = "sample", cols = c( "tomato"), pt.size = 0.5)

DimPlot(srt[,srt$sample == "yg1"], reduction = "umap", group.by = "sample", cols = c( "green3"), pt.size = 0.5)

DimPlot(srt[,srt$sample == "yg3"], reduction = "umap", group.by = "sample", cols = c( "dodgerblue"), pt.size = 0.5)

DimPlot(srt[,srt$sample == "yg5"], reduction = "umap", group.by = "sample", cols = c( "purple"), pt.size = 0.5)

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
##  [3] Biobase_2.54.0              GenomicRanges_1.46.1       
##  [5] GenomeInfoDb_1.30.1         IRanges_2.28.0             
##  [7] S4Vectors_0.32.4            BiocGenerics_0.40.0        
##  [9] MatrixGenerics_1.6.0        matrixStats_0.62.0         
## [11] ggplot2_3.3.6               SeuratObject_4.1.3         
## [13] Seurat_4.1.0               
## 
## loaded via a namespace (and not attached):
##   [1] plyr_1.8.7             igraph_1.3.5           lazyeval_0.2.2        
##   [4] sp_1.5-0               splines_4.1.2          listenv_0.8.0         
##   [7] scattermore_0.8        digest_0.6.30          htmltools_0.5.3       
##  [10] fansi_1.0.3            magrittr_2.0.3         tensor_1.5            
##  [13] cluster_2.1.4          ROCR_1.0-11            globals_0.16.1        
##  [16] spatstat.sparse_3.0-0  colorspace_2.0-3       ggrepel_0.9.1         
##  [19] xfun_0.34              dplyr_1.0.10           RCurl_1.98-1.9        
##  [22] jsonlite_1.8.3         progressr_0.11.0       spatstat.data_3.0-0   
##  [25] survival_3.4-0         zoo_1.8-11             glue_1.6.2            
##  [28] polyclip_1.10-4        gtable_0.3.1           zlibbioc_1.40.0       
##  [31] XVector_0.34.0         leiden_0.4.3           DelayedArray_0.20.0   
##  [34] future.apply_1.9.1     abind_1.4-5            scales_1.2.1          
##  [37] DBI_1.1.3              spatstat.random_3.0-1  miniUI_0.1.1.1        
##  [40] Rcpp_1.0.9             viridisLite_0.4.1      xtable_1.8-4          
##  [43] reticulate_1.26        spatstat.core_2.4-4    htmlwidgets_1.5.4     
##  [46] httr_1.4.4             RColorBrewer_1.1-3     ellipsis_0.3.2        
##  [49] ica_1.0-3              farver_2.1.1           pkgconfig_2.0.3       
##  [52] sass_0.4.2             uwot_0.1.14            deldir_1.0-6          
##  [55] utf8_1.2.2             labeling_0.4.2         tidyselect_1.2.0      
##  [58] rlang_1.0.6            reshape2_1.4.4         later_1.3.0           
##  [61] munsell_0.5.0          tools_4.1.2            cachem_1.0.6          
##  [64] cli_3.4.1              generics_0.1.3         ggridges_0.5.4        
##  [67] evaluate_0.17          stringr_1.4.1          fastmap_1.1.0         
##  [70] yaml_2.3.6             goftest_1.2-3          knitr_1.40            
##  [73] fitdistrplus_1.1-8     purrr_0.3.5            RANN_2.6.1            
##  [76] pbapply_1.5-0          future_1.28.0          nlme_3.1-160          
##  [79] mime_0.12              compiler_4.1.2         rstudioapi_0.14       
##  [82] plotly_4.10.0          png_0.1-7              spatstat.utils_3.0-1  
##  [85] tibble_3.1.8           bslib_0.4.0            stringi_1.7.8         
##  [88] highr_0.9              lattice_0.20-45        Matrix_1.5-1          
##  [91] vctrs_0.5.0            pillar_1.8.1           lifecycle_1.0.3       
##  [94] spatstat.geom_3.0-3    lmtest_0.9-40          jquerylib_0.1.4       
##  [97] RcppAnnoy_0.0.19       data.table_1.14.4      cowplot_1.1.1         
## [100] bitops_1.0-7           irlba_2.3.5.1          httpuv_1.6.6          
## [103] patchwork_1.1.2        R6_2.5.1               promises_1.2.0.1      
## [106] KernSmooth_2.23-20     gridExtra_2.3          parallelly_1.32.1     
## [109] codetools_0.2-18       MASS_7.3-58.1          assertthat_0.2.1      
## [112] withr_2.5.0            sctransform_0.3.5      GenomeInfoDbData_1.2.7
## [115] mgcv_1.8-41            parallel_4.1.2         grid_4.1.2            
## [118] rpart_4.1.19           tidyr_1.2.1            rmarkdown_2.17        
## [121] Rtsne_0.16             shiny_1.7.2